%
% Try Relevant Subtask Learning based classification with
% variational Bayes optimization.
% Note: Ignores images 10-11 as some were missing
%

n_tasks = 9;
pca_dimensionality = 30;

% compute mean vector and PCA projection from all data
x = [];
for img=1:n_tasks;
  x=[x;males{img}];
end;
for img=1:n_tasks;
  x=[x;females{img}];
end;
datamean=mean(x,1);
datacov=cov(x);
[pcaproj,D]=eig(datacov);
d=diag(D); [s,I]=sort(-d);
pcaproj=pcaproj(:,I(1:pca_dimensionality));

%[pcaproj,data] = princomp(x);
%pcaproj=pcaproj(:,1:pca_dimensionality);


% transform data into format required by RSL code: [features 1 classlabel] for each image
tasksmales = cell(n_tasks,1);
tasksfemales = cell(n_tasks,1);
for i=1:n_tasks,
  % Substract overall data mean,
  % apply PCA,
  % append a constant feature (value 1) to the end of the data,
  % and append the class (here -1 males, +1 females)
  tempdata1=males{i};
  tempdata1=[(tempdata1-repmat(datamean,[size(tempdata1,1) 1]))*pcaproj ones(size(tempdata1,1),1) -1*ones(size(tempdata1,1),1)];
  tasksmales{i} = tempdata1;
  tempdata2=females{i};
  tempdata2=[(tempdata2-repmat(datamean,[size(tempdata1,1) 1]))*pcaproj ones(size(tempdata2,1),1) 1*ones(size(tempdata2,1),1)];
  tasksfemales{i} = tempdata2;
end;

% create a stratified random permutation of the 40 persons;
% leave 5 males and 5 females in training sets, rest in test sets
n_trainpersons=7;
traintasks=cell(n_tasks,1);
testtasks=cell(n_tasks,1);
for i=1:n_tasks,
  ord1 = randperm(20);
  ord2 = randperm(20);
  traintasks{i}=[tasksmales{i}(ord1(1:n_trainpersons),:); tasksfemales{i}(ord2(1:n_trainpersons),:)];
  testtasks{i}=[tasksmales{i}(ord1((n_trainpersons+1):end),:); tasksfemales{i}(ord2((n_trainpersons+1):end),:)];
end;


outputfilter=cell(3,n_tasks);

 
%-----------------------------------------------------
% Initialize, run and evaluate VB single-task model
%-----------------------------------------------------

stlqualities = zeros(n_tasks,1);
for task=1:n_tasks,
  % train using only the task-of-interest
  % 0.05 is the prior variance of parameters
  [posterior,prior] = initialize_rsl({traintasks{task}}, 0.05); 
  % 500 is number of iterations
  posterior2 = optimize_rsl(posterior, prior, {traintasks{task}}, 500); 
  classprobabilities = predict_rsl(posterior2, testtasks{task});
  classifications = -1 + 2*(classprobabilities > 0.5);
  classaccuracy = sum(classifications == testtasks{task}(:,end))/size(testtasks{task},1);
  stlqualities(task) = classaccuracy;

  outputfilter{1,task} = posterior2{2}(1:pca_dimensionality);
end;


%-----------------------------------------------------
% Initialize, run and evaluate pooled VB single-task model
%-----------------------------------------------------

poolqualities = zeros(n_tasks,1);
for task=1:n_tasks,
  % create a pooled data set
  tempdata = [];
  for t=1:n_tasks,
    tempdata=[tempdata;traintasks{t}];
  end;

  % train using the pooled data set
  % 0.05 is the prior variance of parameters
  [posterior,prior] = initialize_rsl({tempdata}, 0.05); 
  % 500 is number of iterations
  posterior2 = optimize_rsl(posterior, prior, {tempdata}, 500);  
  classprobabilities = predict_rsl(posterior2, testtasks{task});
  classifications = -1 + 2*(classprobabilities > 0.5);
  classaccuracy = sum(classifications == testtasks{task}(:,end))/size(testtasks{task},1);
  poolqualities(task) = classaccuracy;

  outputfilter{2,task} = posterior2{2}(1:pca_dimensionality);
end;



%-----------------------------------------------------
% Initialize, run and evaluate VB RSL model
%-----------------------------------------------------

rslqualities = zeros(n_tasks,1);
for task=1:n_tasks,

  % switch the task-of-interest to be the first one...
  temptrain = traintasks;
  temptrain{1} = traintasks{task};
  temptrain{task} = traintasks{1};

  % train RSL with the switched tasks
  % 0.05 is the prior variance of parameters
  [posterior,prior] = initialize_rsl(temptrain, 0.05); 
  % 500 is number of iterations
  posterior2 = optimize_rsl(posterior, prior, temptrain, 500); 

  % predict for test data and evaluate
  classprobabilities = predict_rsl(posterior2, testtasks{task});
  classifications = -1 + 2*(classprobabilities > 0.5);
  classaccuracy = sum(classifications == testtasks{task}(:,end))/size(testtasks{task},1);
  rslqualities(task) = classaccuracy;

  outputfilter{3,task} = posterior2{2}(1:pca_dimensionality);
end;


figure;
for method=1:3,
  for task=1:n_tasks,
    subplot(3,n_tasks,(method-1)*n_tasks+task);
    imagesc(reshape(pcaproj*outputfilter{method,task},[16 16]));
    colormap(gray);
  end;
end;
rslqualities
poolqualities
stlqualities

